home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / D-G / GemsII / viewcorr / matrix.c next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  4.6 KB  |  170 lines  |  [TEXT/MPS ]
  1. /*
  2.  * matrix.c - Simple routines for general sized matrices.
  3.  *   
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include "matrix.h"
  9.  
  10.  
  11. double 
  12. InvertMatrix(mat,actual_size)
  13. Matrix mat;            /* Holds the original and inverse */
  14. int actual_size;    /* Actual size of matrix in use, (high_subscript+1)*/
  15. {
  16.     int i,j,k;
  17.                     /* Locations of pivot elements */
  18.     int *pvt_i, *pvt_j;
  19.     double pvt_val;                     /* Value of current pivot element */
  20.     double hold;                        /* Temporary storage */
  21.     double determ;                      /* Determinant */
  22.  
  23.     determ = 1.0;
  24.  
  25.     pvt_i = (int *) malloc(actual_size * sizeof(int));
  26.     pvt_j = (int *) malloc(actual_size * sizeof(int));
  27.  
  28.     for (k = 0; k < actual_size; k++)
  29.     {
  30.         /* Locate k'th pivot element */
  31.         pvt_val = mat[k][k];            /* Initialize for search */
  32.         pvt_i[k] = k;
  33.         pvt_j[k] = k;
  34.         for (i = k; i < actual_size; i++)
  35.           for (j = k; j < actual_size; j++)
  36.             if (fabs(mat[i][j]) > fabs(pvt_val))
  37.             {
  38.                 pvt_i[k] = i;
  39.                 pvt_j[k] = j;
  40.                 pvt_val = mat[i][j];
  41.             }
  42.         /* Product of pivots, gives determinant when finished */
  43.         determ *= pvt_val;
  44.         if (determ == 0.0) {    
  45.          /* Matrix is singular (zero determinant). */
  46.         free(pvt_i);
  47.         free(pvt_j);
  48.             return (0.0);              
  49.     }
  50.  
  51.         /* "Interchange" rows (with sign change stuff) */
  52.         i = pvt_i[k];
  53.         if (i != k)                     /* If rows are different */
  54.           for (j = 0; j < actual_size; j++)
  55.           {
  56.             hold = -mat[k][j];
  57.             mat[k][j] = mat[i][j];
  58.             mat[i][j] = hold;
  59.           }
  60.  
  61.         /* "Interchange" columns */
  62.         j = pvt_j[k];
  63.         if (j != k)                     /* If columns are different */
  64.           for (i = 0; i < actual_size; i++)
  65.           {
  66.             hold = -mat[i][k];
  67.             mat[i][k] = mat[i][j];
  68.             mat[i][j] = hold;
  69.           }
  70.         /* Divide column by minus pivot value */
  71.         for (i = 0; i < actual_size; i++)
  72.           if (i != k)                   /* Don't touch the pivot entry */
  73.             mat[i][k] /= ( -pvt_val) ;  /* (Tricky C syntax for division) */
  74.  
  75.         /* Reduce the matrix */
  76.         for (i = 0; i < actual_size; i++)
  77.         {
  78.             hold = mat[i][k];
  79.             for (j = 0; j < actual_size; j++)
  80.               if ( i != k && j != k )   /* Don't touch pivot. */
  81.                 mat[i][j] += hold * mat[k][j];
  82.         }
  83.  
  84.         /* Divide row by pivot */
  85.         for (j = 0; j < actual_size; j++)
  86.           if (j != k)                   /* Don't touch the pivot! */
  87.             mat[k][j] /= pvt_val;
  88.  
  89.         /* Replace pivot by reciprocal (at last we can touch it). */
  90.         mat[k][k] = 1.0/pvt_val;
  91.     }
  92.  
  93.     /* That was most of the work, one final pass of row/column interchange */
  94.     /* to finish */
  95.     for (k = actual_size-2; k >= 0; k--)  /* Don't need to work with 1 by 1 */
  96.                                         /* corner */
  97.     {
  98.         i = pvt_j[k];         /* Rows to swap correspond to pivot COLUMN */
  99.         if (i != k)                     /* If rows are different */
  100.           for(j = 0; j < actual_size; j++)
  101.           {
  102.             hold = mat[k][j];
  103.             mat[k][j] = -mat[i][j];
  104.             mat[i][j] = hold;
  105.           }
  106.  
  107.         j = pvt_i[k];           /* Columns to swap correspond to pivot ROW */
  108.         if (j != k)                     /* If columns are different */
  109.           for (i = 0; i < actual_size; i++)
  110.           {
  111.             hold = mat[i][k];
  112.             mat[i][k] = -mat[i][j];
  113.             mat[i][j] = hold;
  114.           }
  115.     }
  116.  
  117.     free(pvt_i);
  118.     free(pvt_j);
  119.     return(determ);
  120. }
  121.  
  122. Matrix
  123. NewMatrix(cols, rows)
  124. int cols,rows;
  125. {
  126.     int i;
  127.     Matrix newM;
  128.     newM = (double **) malloc(rows * sizeof(double *));
  129.     for(i = 0; i < rows; i++)
  130.     newM[i] = (double *) malloc(cols * sizeof(double));
  131.     return newM;
  132. }
  133.  
  134. FreeMatrix(mat, rows)
  135. Matrix mat;
  136. int rows;
  137. {
  138.     int i;
  139.     for(i = 0; i < rows; i++)
  140.     free(mat[i]);
  141.     free(mat);
  142. }
  143.  
  144. TransposeMatrix(inM, outM, cols, rows)
  145. Matrix inM, outM;
  146. int cols,rows;
  147. {
  148.     int tempI, tempJ;
  149.     for(tempI=0; tempI < rows; tempI++)
  150.     for(tempJ=0; tempJ < cols; tempJ++)
  151.         outM[tempI][tempJ] = inM[tempJ][tempI];    
  152. }
  153.  
  154. MultMatrix(firstM, secondM, outM, firstrows, cols, secondcols)
  155. Matrix firstM, secondM, outM;
  156. int firstrows, cols, secondcols;
  157. {
  158.     int i,j,k;
  159.     double sum;
  160.  
  161.     for(i=0; i < secondcols; i++)
  162.     for(j=0; j < firstrows; j++)
  163.     {
  164.         sum = 0.0;
  165.         for(k=0; k < cols; k++)
  166.         sum += firstM[j][k] * secondM[k][i];
  167.         outM[j][i] = sum;
  168.     }
  169. }
  170.